% Find SELECT for changes to sample, specification, solver, etc

clc;
clear;
format long

% Import data file, SELECT sample, and label variables:
data=importdata('Estimation_data_all.txt');                 % N=2682: all completed main-draw matches with non-missing environmental conditions

% SELECT ability_metric 1=pre-match odds, 2.5=granular rank, or 4=z-scores
ability_metric=1;                    % E.g., Table 5, column 1; Table A.4, column 1
%ability_metric=2.5;                 % E.g., Table 5, column 2; Table A.4, column 2
%ability_metric=4;                   % E.g., Table 5, column 3; Table A.4, column 3
% SELECT prize spread 1=ones, 2=cash prize spread, 3=points prize spread
prize_metric=2;
% SELECT finite difference for covariance
h=1e-02;


augmented_contest_prize=0;
smallDiff=data(:,14);                             % N=700 if further restrict to symmetric matches (smallDiff is indicator for winProb_diff<=30)
closeFirstSet=(data(:,1)+data(:,2)>=10);          % Alternative close first set to address Malueg and Yates concern (N=320, symmetric w/ ability_metric=1): 6-4, 4-6, 5-5, 6-5 or 5-6 
ageSample=data(:,23);                                                                                % age is available for both players
%ageSample=(data(:,23)==1&data(:,26)<=3.6);                                                          % age restriction 1: age difference below median of 3.6y (besides age is available for both players)    
%ageSample=(data(:,23)==1&data(:,26)>3.6);                                                           % age restriction 2: age difference above median of 3.6y (besides age is available for both players)
%ageSample=(data(:,23)==1&data(:,24)<=24.6&data(:,25)<=24.6);                                        % age restriction 3: two below-median-age players (besides age is available for both players)
%ageSample=(data(:,23)==1&data(:,24)>24.6&data(:,25)>24.6);                                          % age restriction 4: two over-median-age players (besides age is available for both players)
%ageSample=(data(:,23)==1&((data(:,24)<=24.6&data(:,25)>24.6)|(data(:,24)>24.6 &data(:,25)<=24.6))); % age restriction 5: one below-median-age player and one over-median-age player (besides age is available for both players)
dvRain=data(:,31);                                % record of rain is available for match
notHugePrize=(data(:,27)<=25.896);                % matches with prize no greater than the 75th percentile over 2811 main-draw matches
%prestige=(data(:,33)==1|data(:,34)==1);               % Australian Open, China Open
%prestige=(data(:,33)==1|data(:,34)==1|data(:,35)==1); % Australian Open, China Open, Wuhan Open
prestige=(data(:,33)==0);                              % Only matches in China
round=data(:,37);
%roofClosed=data(:,38);                           % 13 completed matches played under closed roof due to heat. Only 3 after 2008, when we observe betting odds.
%roofClosed=data(:,39);                           % 53 completed matches played under closed roof due to rain. 45 after 2008, when we observe betting odds.
roofClosed=(data(:,38)==1|data(:,39)==1);

if ability_metric==1
    dvBet=data(:,10);                              % 1. pre-match odds are available for both players
    data_sample=data(dvBet==1,:);                  % N=2157: all among the 2682 matches for which pre-match odds are available (both players)
    %augmented_contest_prize=1;                    % (Robustness Table 9) ACTIVATE to run augmented contest model variant
    %data_sample=data(dvBet==1&smallDiff==1,:);    % (Robustness: Table 8, column 5) SELECT symmetric (w/ ability_metric=1). N=700: the 700 symmetric matches; dvBet==1 is redundant since pre-match odds are available (both players)
    %data_sample=data(dvBet==1&smallDiff==1&closeFirstSet==1,:);  % (Robustness: Footnote 36 of the article) SELECT symmetric & close first set
    %data_sample=data(dvBet==1&round>=8,:);        % (Robustness: Table 8, column 2) SELECT sample without quarterfinal (R4), semifinal (R2), and final (R1)
    %data_sample=data(dvBet==1&ageSample==1,:);    % SELECT if player ages or age difference is used (missing for one match) 
    %data_sample=data(dvBet==1&notHugePrize==1,:); % (Robustness: Table 8, column 3) SELECT if drop contests with huge prizes (as defined above)
    %data_sample=data(dvBet==1&prestige==1,:);     % (Robustness: Table 8, column 4) SELECT if drop contests with less prestige
    %data_sample=data(dvBet==1&roofClosed==0,:);   % (Robustness: Table A.4, column 1) SELECT if drop matches plausibly played under closed roof due to heat
elseif ability_metric==2.5
    dvRank=data(:,15);                             % 2.5. (granular) WTA ranks are available for both players
    data_sample=data(dvRank==1,:);                 % N=2670: all among the 2682 matches for which WTA ranks are available (both players)
    %data_sample=data(dvRank==1&ageSample==1,:);   % SELECT if player ages or age difference is used (missing for one match) 
    %data_sample=data(dvRank==1&roofClosed==0,:);  % (Robustness: Table A.4, column 2) SELECT if drop matches plausibly played under closed roof due to heat
elseif ability_metric==4
    dvPts=data(:,18);                              % 4. normalized ranking points (z-scores) are available for both players
    data_sample=data(dvPts==1,:);                  % N=2670: all among the 2682 matches for which WTA ranking points are available (both players)--this happens to be the same sample as in 2
    %data_sample=data(dvPts==1&ageSample==1,:);    % SELECT if player ages or age difference is used (missing for one match) 
    %data_sample=data(dvPts==1&roofClosed==0,:);   % (Robustness: Table A.4, column 3) SELECT if drop matches plausibly played under closed roof due to heat
end

clear data dvBet dvRank dvPts smallDiff round closeFirstSet ageSample dvRain notHugePrize prestige roofClosed;

N=size(data_sample,1);

w1=data_sample(:,1);           % set 1 games won by match winner
l1=data_sample(:,2);           % set 1 games won by match loser
w2=data_sample(:,3);           % set 2 games won by match winner
l2=data_sample(:,4);           % set 2 games won by match loser
w3=data_sample(:,5);           % set 3 games won by match winner (if applicable)
l3=data_sample(:,6);           % set 3 games won by match loser (if applicable)
%wSets=data_sample(:,7);       % number of sets won by match winner
%lSets=data_sample(:,8);       % number of sets won by match loser
%tSets=data_sample(:,9);       % total number of sets
wWinProb=data_sample(:,11);    % (if pre-match odds are available) ex-ante winning probability of match winner
lWinProb=data_sample(:,12);    % (if pre-match odds are available) ex-ante winning probability of match loser
winProbDiff=data_sample(:,13); % abs(wWinProb-lWinProb)
wRank=data_sample(:,16);       % (if both players are ranked) WTA rank of match winner
lRank=data_sample(:,17);       % (if both players are ranked) WTA rank of match loser
%wPts=data_sample(:,19);       % (if both players have points) WTA ranking points of match winner
%lPts=data_sample(:,20);       % (if both players have points) WTA ranking points of match loser
wPtsN=data_sample(:,21);       % (if both players have points) WTA z-score of match winner
lPtsN=data_sample(:,22);       % (if both players have points) WTA z-score of match loser
%wAge=data_sample(:,24);       % (if age is available) age of match winner
%lAge=data_sample(:,25);       % (if age is available) age of match loser
%ageDiff=data_sample(:,26);    % (if age is available) abs(wAge-lAge)
prizeCashK=data_sample(:,27);  % win-lose cash prize spread of match (in 1000 USD)
prizePts=data_sample(:,28);    % win-lose WTA ranking points prize spread of match  
tpORtp3h=data_sample(:,29);    % temperature during match
Cpm25=data_sample(:,30);       % pollution during match
raining=data_sample(:,32);	   % (if record of rain is available) dummy for any rain during match
%dvAussieO=data_sample(:,33);  % Australian Open match
%dvChinaO=data_sample(:,34);   % China Open match
%dvWuhanO=data_sample(:,35);   % Wuhan Open match
%dvOtherC=data_sample(:,36);   % Other Chinese WTA series match (other than in Beijing and in Wuhan)
round=data_sample(:,37);       % round of WTA series
cityID=data_sample(:,40);      % cityID of WTA series
yearSeries=data_sample(:,41);  % year of WTA series

clear data_sample;
       

% We want to label players according to low (marginal) cost (of effort) "l" and high cost "h" (not winner "w" and loser "l")
% Cost will be a function of ability, so index match winner according to whether she was the low-cost (equivalently, more-able) player

if ability_metric==1
    LowCost_WinsMatch=(wWinProb>=lWinProb);                                                             % (1596 out of 2157 matches = 74%) 1. Ability proxy is pre-match winning probability (betting odds)
    LowCost_Ability =(.01*wWinProb).*LowCost_WinsMatch + (.01*lWinProb).*(ones(N,1)-LowCost_WinsMatch); % max(wWinProb,lWinProb) and min(wWinProb,lWinProb), respectively, would work just as well
    HighCost_Ability=(.01*lWinProb).*LowCost_WinsMatch + (.01*wWinProb).*(ones(N,1)-LowCost_WinsMatch); % (used only for observations where pre-match odds are available for both players)
elseif ability_metric==2.5
    LowCost_WinsMatch=(wRank<=lRank);                                                                   % (1871 out of 2670 matches = 70%) 2.5. (granular) Ability proxy is WTA rank
    LowCost_Ability =wRank.*LowCost_WinsMatch    + lRank.*(ones(N,1)-LowCost_WinsMatch);                % min(wRank,lRank) and max(wRank,lRank), respectively, would work just as well
    HighCost_Ability=lRank.*LowCost_WinsMatch    + wRank.*(ones(N,1)-LowCost_WinsMatch);                % (used only for observations where WTA ranks are available for both players)
elseif ability_metric==4
    LowCost_WinsMatch=(wPtsN>=lPtsN);                                                                   % (1870 out of 2670 matches = 70%) 4. Ability proxy is WTA z-score (normalized ranking points)
    LowCost_Ability =wPtsN.*LowCost_WinsMatch    + lPtsN.*(ones(N,1)-LowCost_WinsMatch);                % max(wPtsN,lPtsN) and min(wPtsN,lPtsN), respectively, would work just as well
    HighCost_Ability=lPtsN.*LowCost_WinsMatch    + wPtsN.*(ones(N,1)-LowCost_WinsMatch);                % (used only for observations where WTA ranking points are available for both players)   
end


% Prepare indicators for match outcomes
LowCost_Set1Games =w1.*(LowCost_WinsMatch==1)+l1.*(LowCost_WinsMatch==0);
HighCost_Set1Games=l1.*(LowCost_WinsMatch==1)+w1.*(LowCost_WinsMatch==0);
LowCost_Set2Games =w2.*(LowCost_WinsMatch==1)+l2.*(LowCost_WinsMatch==0);
HighCost_Set2Games=l2.*(LowCost_WinsMatch==1)+w2.*(LowCost_WinsMatch==0);
LowCost_Set3Games =w3.*(LowCost_WinsMatch==1)+l3.*(LowCost_WinsMatch==0);
HighCost_Set3Games=l3.*(LowCost_WinsMatch==1)+w3.*(LowCost_WinsMatch==0);
outcome_ll= ( (LowCost_Set1Games>HighCost_Set1Games) & (LowCost_Set2Games>HighCost_Set2Games) );
outcome_lhl=( (LowCost_Set1Games>HighCost_Set1Games) & (LowCost_Set2Games<HighCost_Set2Games) & (LowCost_Set3Games>HighCost_Set3Games) );
outcome_lhh=( (LowCost_Set1Games>HighCost_Set1Games) & (LowCost_Set2Games<HighCost_Set2Games) & (LowCost_Set3Games<HighCost_Set3Games) );
outcome_hll=( (LowCost_Set1Games<HighCost_Set1Games) & (LowCost_Set2Games>HighCost_Set2Games) & (LowCost_Set3Games>HighCost_Set3Games) );
outcome_hlh=( (LowCost_Set1Games<HighCost_Set1Games) & (LowCost_Set2Games>HighCost_Set2Games) & (LowCost_Set3Games<HighCost_Set3Games) );
outcome_hh= ( (LowCost_Set1Games<HighCost_Set1Games) & (LowCost_Set2Games<HighCost_Set2Games) );
outcome=[outcome_ll outcome_lhl outcome_lhh outcome_hll outcome_hlh outcome_hh];
clear LowCost_Set1Games HighCost_Set1Games LowCost_Set2Games HighCost_Set2Games LowCost_Set3Games HighCost_Set3Games;
clear outcome_ll outcome_lhl outcome_lhh outcome_hll outcome_hlh outcome_hh;


% w/l in original variable names is confusing, since l denotes loser, not low-cost player; thus cleanse the variables that will not be used further
clear wRank lRank wPts lPts wPtsN lPtsN l1 l2 l3 w1 w2 w3;
clear wAge lAge ageDiff;

if prize_metric==1
    V_l=ones(N,1);
    V_h=V_l;
elseif prize_metric==2
    V_l=prizeCashK;
    V_h=V_l;
elseif prize_metric==3
    V_l=prizePts;
    V_h=V_l;
end
clear prizeCash prizePts;

% Parameters that are fixed
cutoff_tp=27;                           % Coefficient can apply on bin or linearly in the difference or log difference beyond threshold
cutoff_pm=150/100;
delta_0=0;

% Parameters that are estimated: Initial values over which search takes place.
k_ini       =1;
if ability_metric==1
    theta_ini=0;            % 1. Ability proxy is pre-match winning probability (betting odds)
elseif ability_metric==2.5
    theta_ini=0;            % 2.5. Ability proxy is WTA rank (granular)
elseif ability_metric==4
    theta_ini=0;            % 4. Ability proxy is WTA z-score (normalized ranking points)
end
delta_tp_ini=0; 
delta_pm_ini=0;

parameter_fixed=[ cutoff_tp;            % Temperature (tpORtp3h) cutoff
                  cutoff_pm;            % PM2.5 (Cpm25) cutoff
                  delta_0 ];            % Delta intercept
parameter_ini=[ k_ini;                  % technology parameter (randomness with which cost differences translate into winning probabilities)
                theta_ini;              % cost parameter
                delta_tp_ini;           % temperature disutility parameter (variation in sample is 0 = 25 C to 19 = 44 C)
                delta_pm_ini ];         % PM2.5 disutility parameter (variation in sample is 0 = 100 ug/m3 to 4.2 = 520 ug/m3)
clear k_ini theta_ini delta_tp_ini delta_pm_ini;
        
disp('The mean (negative of the) log likelihood contribution across matches at the initial parameter values is ')
likelihood_obj(parameter_ini,ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed)


% Optimization
% knitro with constraints
options = optimset('algorithm','interior-point','display','iter','maxiter',100000,'maxfunevals',100000,'TolX',1e-20,'TolFun',1e-20,'TolCon',1e-6,'Hessian','bfgs');
A=[]; b=[];
lb=[0; 0;   -Inf; -Inf];
ub=[1; Inf; Inf;  Inf];
tic
%[parameter,fval,exitflag,output]=knitromatlab(@(parameter) likelihood_obj(parameter,ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed),parameter_ini,A,b,[],[],lb,ub,[],[],options,[])
[parameter,fval,exitflag,output]=knitromatlab(@(parameter) likelihood_obj(parameter,ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed),parameter_ini,[],[],[],[],lb,ub,[],[],options,[])
%[parameter,fval,exitflag,output]=knitromatlab(@(parameter) likelihood_obj(parameter,ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed),parameter_ini,[],[],[],[],[],[],[],[],options,[])
toc

if ability_metric==4                                            % 4. Ability proxy is WTA z-score (normalized ranking points)
    parameter_ini=[parameter(1); parameter(2); .1; 1.1];        % One more optimization, using estimates of environmental parameters based on betting odds as initial values
    [parameter,fval,exitflag,output]=knitromatlab(@(parameter) likelihood_obj(parameter,ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed),parameter_ini,[],[],[],[],lb,ub,[],[],options,[])
    % Confirm with fminsearch
    options=optimset('display','iter','maxiter',100000,'maxfunevals',100000,'TolX',1e-15,'TolFun',1e-10);
    [parameter,fval,exitflag,output]=fminsearch(@(parameter) likelihood_obj(parameter,ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed), parameter_ini, options)
end


% Robustness to optimizer (see Table 5 caption): fminsearch
%{
options=optimset('display','iter','maxiter',100000,'maxfunevals',100000,'TolX',1e-15,'TolFun',1e-10);
tic
[parameter,fval,exitflag,output]=fminsearch(@(parameter) likelihood_obj(parameter,ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed), parameter_ini, options)
toc
%}



% (Robustness Table 9) This augmented contest prize model variant was written for ability_metric 1's estimation sample defined by dvBet==1 (otherwise need to adjust code)
if (ability_metric==1) && (augmented_contest_prize==1)
    
    % round 1 (final), 2, 4, 8, ..., 64 -> roundSeq 1 (final), 2, 3, 4, ..., 7
    roundSeq = ones(N,1).*(round==1) + repmat(2,N,1).*(round==2) + repmat(3,N,1).*(round==4) + repmat(4,N,1).*(round==8) + repmat(5,N,1).*(round==16) + repmat(6,N,1).*(round==32) + repmat(7,N,1).*(round==64);

    % Takes estimated k from previous optimization, as well as symmetric match prize money, and computes contestprize which accounts for continuation value from subsequent rounds of matches (against symmetric opponent, mild environment)
    prizeMoney=V_l;                 % V_l, to contain the contest prize, will be computed from a city-year series' match-prize-money vector (invariant within round) for every new parameter estimate (specifically k) 

    parameter_000=parameter;
    flag=0;
    count=1;        % One optimization above -> Augmented contest prize -> Optimization
    toler=1e-4;

    while flag==0
        count=count+1

        clear V_l V_h;
        k=parameter(1);

        V_l=zeros(N,1);
        for c=1:7
        for y=2016:-1:2008          % We are missing odds prior to 2008, so no need to go back further    
        exppayoff_nextR=0; 
        for r=1:7                   % For each city-year series, start with the final and work recursively
            prizemoney_cyr=mean(prizeMoney((cityID==c)&(yearSeries==y)&(roundSeq==r)));
                    % The following city-year-round triples had matches held but all matches for that triple are missing in the data (ignores minor year-on-year prize money adjustments unless noted otherwise) 
                    if     (c==2)&&(y==2008)&&(r==1)       %use y+1=2009 * 0.11 factor
                        prizemoney_cyr=mean(prizeMoney((cityID==c)&(yearSeries==y+1)&(roundSeq==r))) * 0.11;
                    elseif (c==2)&&(y==2011)&&(r==4)       %use y-1=2010 (and similarly below)
                        prizemoney_cyr=mean(prizeMoney((cityID==c)&(yearSeries==y-1)&(roundSeq==r)));
                    elseif (c==2)&&(y==2011)&&(r==5)
                        prizemoney_cyr=mean(prizeMoney((cityID==c)&(yearSeries==y-1)&(roundSeq==r)));
                    elseif (c==3)&&(y==2015)&&(r==1)
                        prizemoney_cyr=mean(prizeMoney((cityID==c)&(yearSeries==y-1)&(roundSeq==r)));
                    elseif (c==4)&&(y==2016)&&(r==1)
                        prizemoney_cyr=mean(prizeMoney((cityID==c)&(yearSeries==y-1)&(roundSeq==r)));
                    end
            contestprize=prizemoney_cyr+exppayoff_nextR;        % E.g., in round 2 (semifinal), this is the prize money in round 2 and the expected payoff in round 1 (final), both of which the player enjoys by winning in round 2
            if ~isnan(contestprize)
                V_l((cityID==c)&(yearSeries==y)&(roundSeq==r))=contestprize;
            end
            exppayoff_nextR = 0.5*contestprize - exp_effort_cost(contestprize,k); % E.g. contd., Besides the semifinal's prize money, effort in the semifinal responds to continuation value (in the final) from a win (four lines above: contestprize=prizemoney_cyr+exppayoff_nextR)
                                                                                  % To be used when obtaining expected effort in the quarterfinal: the expected payoff in the semifinal is the contest prize with a 50% likelihood of success less expected effort cost in the semifinal
        end
        end
        end
        %scatter(prizeMoney,V_l);
        V_h=V_l;

        [parameter,fval,exitflag,output]=knitromatlab(@(parameter) likelihood_obj(parameter,ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed),parameter_ini,[],[],[],[],lb,ub,[],[],options,[])
        max(abs(parameter-parameter_000))

        if max(abs(parameter-parameter_000))<toler
            flag=1;                 % Note that this is not a fixed point in V: at each iteration we start from the prize money (data), add continuation value from the final match (using estimate k and player beliefs on symmetry/environment) recursively to form the contest prize, and re-estimate the model 
        end
        parameter_000=parameter;
    end
    clear parameter_000 exppayoff_nextR prizemoney_cyr contestprize

    % Figure A.7
    median(V_l)
    median(prizeMoney)
    V_l_final             =V_l(roundSeq==1);
    V_l_semiquarter       =V_l(roundSeq==2|roundSeq==3);
    V_l_early             =V_l(roundSeq>3);
    prizeMoney_final      =prizeMoney(roundSeq==1);
    prizeMoney_semiquarter=prizeMoney(roundSeq==2|roundSeq==3);
    prizeMoney_early      =prizeMoney(roundSeq>3);
    figure;
    subplot(1,1,1);
    scatter(prizeMoney_early,V_l_early,'red');
    hold on
    scatter(prizeMoney_semiquarter,V_l_semiquarter,'blue');
    hold on
    scatter(prizeMoney_final,V_l_final,'black');
    hold off
    title('Augmenting the contest prize money with continuation value in a series');
    ylabel('Match prize money plus continuation value ($ \times 1000)');
    xlabel('Match prize money ($ \times 1000)');
    legend('Round 64 to round 8 contest','Quarterfinal and semifinal (rounds 4 and 2)','Series final (contest prize = match prize money)');

    figure;
    subplot(1,1,1);
    scatter(prizeMoney_early,V_l_early,'red');
    title('Augmenting the contest prize money with continuation value in a series');
    ylabel('Match prize money plus continuation value ($ \times 1000)');
    xlabel('Match prize money ($ \times 1000)');
    xlim([0 300]);

end
      


% Covariance
prob_i=likelihood_contrib(parameter,ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);  % prob_i (or f) = likelihood_contrib

f_1plush=likelihood_contrib(parameter+[h;0;0;0] ,ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
f_1lessh=likelihood_contrib(parameter+[-h;0;0;0],ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
f_2plush=likelihood_contrib(parameter+[0;h;0;0] ,ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
f_2lessh=likelihood_contrib(parameter+[0;-h;0;0],ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
f_3plush=likelihood_contrib(parameter+[0;0;h;0] ,ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
f_3lessh=likelihood_contrib(parameter+[0;0;-h;0],ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
f_4plush=likelihood_contrib(parameter+[0;0;0;h] ,ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
f_4lessh=likelihood_contrib(parameter+[0;0;0;-h],ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);

% Compute d prob_i / d phi
D1=(f_1plush-f_1lessh)/(2*h);
D2=(f_2plush-f_2lessh)/(2*h);
D3=(f_3plush-f_3lessh)/(2*h);
D4=(f_4plush-f_4lessh)/(2*h);

% Compute d2 prob_i / d phi_j_k
D11=(f_1plush-2*prob_i+f_1lessh)/(h.^2);
D22=(f_2plush-2*prob_i+f_2lessh)/(h.^2);
D33=(f_3plush-2*prob_i+f_3lessh)/(h.^2);
D44=(f_4plush-2*prob_i+f_4lessh)/(h.^2);

f_1plush_2plush=likelihood_contrib(parameter+[h;h;0;0],ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
f_1plush_2lessh=likelihood_contrib(parameter+[h;-h;0;0],ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
f_1lessh_2plush=likelihood_contrib(parameter+[-h;h;0;0],ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
f_1lessh_2lessh=likelihood_contrib(parameter+[-h;-h;0;0],ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
f_1plush_3plush=likelihood_contrib(parameter+[h;0;h;0],ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
f_1plush_3lessh=likelihood_contrib(parameter+[h;0;-h;0],ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
f_1lessh_3plush=likelihood_contrib(parameter+[-h;0;h;0],ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
f_1lessh_3lessh=likelihood_contrib(parameter+[-h;0;-h;0],ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
f_1plush_4plush=likelihood_contrib(parameter+[h;0;0;h],ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
f_1plush_4lessh=likelihood_contrib(parameter+[h;0;0;-h],ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
f_1lessh_4plush=likelihood_contrib(parameter+[-h;0;0;h],ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
f_1lessh_4lessh=likelihood_contrib(parameter+[-h;0;0;-h],ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
D12=((f_1plush_2plush-f_1lessh_2plush)/(2*h)-(f_1plush_2lessh-f_1lessh_2lessh)/(2*h))/(2*h);
D13=((f_1plush_3plush-f_1lessh_3plush)/(2*h)-(f_1plush_3lessh-f_1lessh_3lessh)/(2*h))/(2*h);
D14=((f_1plush_4plush-f_1lessh_4plush)/(2*h)-(f_1plush_4lessh-f_1lessh_4lessh)/(2*h))/(2*h);

f_2plush_3plush=likelihood_contrib(parameter+[0;h;h;0],ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
f_2plush_3lessh=likelihood_contrib(parameter+[0;h;-h;0],ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
f_2lessh_3plush=likelihood_contrib(parameter+[0;-h;h;0],ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
f_2lessh_3lessh=likelihood_contrib(parameter+[0;-h;-h;0],ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
f_2plush_4plush=likelihood_contrib(parameter+[0;h;0;h],ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
f_2plush_4lessh=likelihood_contrib(parameter+[0;h;0;-h],ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
f_2lessh_4plush=likelihood_contrib(parameter+[0;-h;0;h],ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
f_2lessh_4lessh=likelihood_contrib(parameter+[0;-h;0;-h],ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
D21=((f_1plush_2plush-f_1plush_2lessh)/(2*h)-(f_1lessh_2plush-f_1lessh_2lessh)/(2*h))/(2*h);
D23=((f_2plush_3plush-f_2lessh_3plush)/(2*h)-(f_2plush_3lessh-f_2lessh_3lessh)/(2*h))/(2*h);
D24=((f_2plush_4plush-f_2lessh_4plush)/(2*h)-(f_2plush_4lessh-f_2lessh_4lessh)/(2*h))/(2*h);

f_3plush_4plush=likelihood_contrib(parameter+[0;0;h;h],ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
f_3plush_4lessh=likelihood_contrib(parameter+[0;0;h;-h],ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
f_3lessh_4plush=likelihood_contrib(parameter+[0;0;-h;h],ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
f_3lessh_4lessh=likelihood_contrib(parameter+[0;0;-h;-h],ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed);
D31=((f_1plush_3plush-f_1plush_3lessh)/(2*h)-(f_1lessh_3plush-f_1lessh_3lessh)/(2*h))/(2*h);
D32=((f_2plush_3plush-f_2plush_3lessh)/(2*h)-(f_2lessh_3plush-f_2lessh_3lessh)/(2*h))/(2*h);
D34=((f_3plush_4plush-f_3lessh_4plush)/(2*h)-(f_3plush_4lessh-f_3lessh_4lessh)/(2*h))/(2*h);

D41=((f_1plush_4plush-f_1plush_4lessh)/(2*h)-(f_1lessh_4plush-f_1lessh_4lessh)/(2*h))/(2*h);
D42=((f_2plush_4plush-f_2plush_4lessh)/(2*h)-(f_2lessh_4plush-f_2lessh_4lessh)/(2*h))/(2*h);
D43=((f_3plush_4plush-f_3plush_4lessh)/(2*h)-(f_3lessh_4plush-f_3lessh_4lessh)/(2*h))/(2*h);

H=zeros(4,4,N);

H(1,1,:) = D11./prob_i - D1.*D1./prob_i./prob_i;
H(1,2,:) = D12./prob_i - D1.*D2./prob_i./prob_i;
H(1,3,:) = D13./prob_i - D1.*D3./prob_i./prob_i;
H(1,4,:) = D14./prob_i - D1.*D4./prob_i./prob_i;

H(2,1,:) = D21./prob_i - D2.*D1./prob_i./prob_i;
H(2,2,:) = D22./prob_i - D2.*D2./prob_i./prob_i;
H(2,3,:) = D23./prob_i - D2.*D3./prob_i./prob_i;
H(2,4,:) = D24./prob_i - D2.*D4./prob_i./prob_i;

H(3,1,:) = D31./prob_i - D3.*D1./prob_i./prob_i;
H(3,2,:) = D32./prob_i - D3.*D2./prob_i./prob_i;
H(3,3,:) = D33./prob_i - D3.*D3./prob_i./prob_i;
H(3,4,:) = D34./prob_i - D3.*D4./prob_i./prob_i;

H(4,1,:) = D41./prob_i - D4.*D1./prob_i./prob_i;
H(4,2,:) = D42./prob_i - D4.*D2./prob_i./prob_i;
H(4,3,:) = D43./prob_i - D4.*D3./prob_i./prob_i;
H(4,4,:) = D44./prob_i - D4.*D4./prob_i./prob_i;

Hess=-mean(H,3);
cov_hess=inv(Hess)/N;
se_hess=diag(cov_hess).^(1/2);
tstat_hess=parameter./se_hess;
format shortG
disp('estimate      std dev (Hessian estimate)     t-stat (Hessian estimate)');
[parameter se_hess tstat_hess]








%{
% Figures 8, 9, and A.5
k        =parameter(1);
theta    =parameter(2);
delta_tp =parameter(3);
delta_pm =parameter(4);
shift_tp =log( max( (32      - cutoff_tp) , 0) + 1);        % Shift 5 = 27 to 32 Celsius, regressor increases by 1.7918 
shift_tp2=log( max( (37      - cutoff_tp) , 0) + 1);        % Shift 10 = 27 to 37 Celsius, regressor increases by 2.3979 
shift_pm =log( max( (200/100 - cutoff_pm) , 0) + 1);        % Shift 50 = 150 to 200 ug/m3, regressor increases by 0.4055
shift_pm2=log( max( (250/100 - cutoff_pm) , 0) + 1);        % Shift 100 = 150 to 250 ug/m3, regressor increases by 0.6932

delta_no_gridpoints=5;
delta_dis=[0; delta_tp*shift_tp; delta_tp*shift_tp2; delta_pm*shift_pm; delta_pm*shift_pm2];

if ability_metric==1              % 1. Ability proxy is pre-match winning probability (betting odds): I want a difference in opponents' winning probability varying over 0, 10, ..., 80. 
    cost_no_gridpoints=17;        % Note that in the data, mean(wWinProb+lWinProb) = 106.79, so construct LowCost_Ability_dis & HighCost_Ability_dis accordingly.
    LowCost_Ability_dis =.01*[106.79/2; 106.79/2+2.5; 106.79/2+5; 106.79/2+7.5; 106.79/2+10; 106.79/2+12.5; 106.79/2+15; 106.79/2+17.5; 106.79/2+20; 106.79/2+22.5; 106.79/2+25; 106.79/2+27.5; 106.79/2+30; 106.79/2+32.5; 106.79/2+35; 106.79/2+37.5; 106.79/2+40];
    HighCost_Ability_dis=.01*[106.79/2; 106.79/2-2.5; 106.79/2-5; 106.79/2-7.5; 106.79/2-10; 106.79/2-12.5; 106.79/2-15; 106.79/2-17.5; 106.79/2-20; 106.79/2-22.5; 106.79/2-25; 106.79/2-27.5; 106.79/2-30; 106.79/2-32.5; 106.79/2-35; 106.79/2-37.5; 106.79/2-40];
    winProbDiff_dis=(1/.01)*(LowCost_Ability_dis-HighCost_Ability_dis);
    c_l_dis=1* LowCost_Ability_dis.^(-1*theta);
    c_h_dis=1*HighCost_Ability_dis.^(-1*theta);    
elseif ability_metric==2.5         % 2.5. Ability proxy is WTA rank: I want a difference in opponents' WTA rank of 1 v 1, 1 v 10, 1 v 100, 1 v 1000. Maybe include 10 v 100 etc.
    %{
    cost_no_gridpoints=12;
    LowCost_Ability_dis =[1; 1; 1;  1;  1;  1;  1;  1;   1;   1;   1;    1];
    HighCost_Ability_dis=[1; 2; 5; 10; 20; 30; 40; 50; 100; 200; 500; 1000];
    %}
    cost_no_gridpoints=1000;
    LowCost_Ability_dis =ones(cost_no_gridpoints,1);
    HighCost_Ability_dis=zeros(cost_no_gridpoints,1);
    for i=1:cost_no_gridpoints
        HighCost_Ability_dis(i)=i;
    end
    %c_l_dis=1*(floor( LowCost_Ability_dis/10)+1).^(theta);
    %c_h_dis=1*(floor(HighCost_Ability_dis/10)+1).^(theta);
    c_l_dis =1* LowCost_Ability_dis.^theta;
    c_h_dis =1*HighCost_Ability_dis.^theta;
end

cl_ch_ratio  =c_l_dis./c_h_dis;
Vin          =15.092;    %Figure 8: Median USD 15,092 across 2811 main-draw matches
%Vin          =8.672;    %Figures 9 and A.5: p25 USD 8,672 across 2811 main-draw matches
momentum_l   =zeros(cost_no_gridpoints,delta_no_gridpoints);
prob3setter  =zeros(cost_no_gridpoints,delta_no_gridpoints);
winProb_l_dis=zeros(cost_no_gridpoints,delta_no_gridpoints);
exppayoff_l  =zeros(cost_no_gridpoints,delta_no_gridpoints);
exppayoff_h  =zeros(cost_no_gridpoints,delta_no_gridpoints);
ones         =ones(delta_no_gridpoints,1);
for i=1:cost_no_gridpoints     
    % Each function call fixes player asymmetry and varies environment within call
    [p_llose3,p_llose2_lwon1,p_llose2_llost1,p_llose1,exppayoff_l_case,exppayoff_h_case] = diagnostics(delta_dis,cl_ch_ratio(i),1,Vin,Vin,k);     % Inputs are delta, c_l, c_h, V_l, V_h, k   
    momentum_l(i,:)     =(p_llose2_llost1-p_llose2_lwon1)';     % This is equal to momentum_h
    prob3setter(i,:)    =( (ones-p_llose1).*p_llose2_lwon1.*(ones-p_llose3) + (ones-p_llose1).*p_llose2_lwon1.*p_llose3 + p_llose1.*(ones-p_llose2_llost1).*(ones-p_llose3) + p_llose1.*(ones-p_llose2_llost1).*p_llose3 )';
    winProb_l_dis(i,:)  =( (ones-p_llose1).*(ones-p_llose2_lwon1) + (ones-p_llose1).*p_llose2_lwon1.*(ones-p_llose3) + p_llose1.*(ones-p_llose2_llost1).*(ones-p_llose3) )';
    exppayoff_l(i,:)    =exppayoff_l_case';
    exppayoff_h(i,:)    =exppayoff_h_case';
end
clear exppayoff_l_case exppayoff_h_case


% Low-cost player's set 2 winning probabilities: State dependence
% SELECT: winProbDiff_dis, HighCost_Ability_dis OR cl_ch_ratio along the horizontal axis
figure;
subplot(1,1,1);
plot(winProbDiff_dis,momentum_l(:,1),'g',winProbDiff_dis,momentum_l(:,2),'--r',winProbDiff_dis,momentum_l(:,3),'--ro',winProbDiff_dis,momentum_l(:,4),':k',winProbDiff_dis,momentum_l(:,5),':ko','LineWidth',3);
%plot(HighCost_Ability_dis,momentum_l(:,1),'g',HighCost_Ability_dis,momentum_l(:,2),'--r',HighCost_Ability_dis,momentum_l(:,3),'--ro',HighCost_Ability_dis,momentum_l(:,4),':k',HighCost_Ability_dis,momentum_l(:,5),':ko','LineWidth',3);
%plot(cl_ch_ratio    ,momentum_l(:,1),'g',cl_ch_ratio    ,momentum_l(:,2),'--r',cl_ch_ratio    ,momentum_l(:,3),'--ro',cl_ch_ratio    ,momentum_l(:,4),':k',cl_ch_ratio    ,momentum_l(:,5),':ko','LineWidth',3);
ylabel('State-induced change in battle-2 success rate');
xlabel('Difference in opponents'' pre-match winning probability (%)');
%xlabel('High-cost player''s rank (competing with player ranked 1)');
%xlabel('Players'' relative marginal costs, c_l/c_h');
legend('27 ^{\circ}C, 150 \mug/m^{3}','32 ^{\circ}C, 150 \mug/m^{3}','37 ^{\circ}C, 150 \mug/m^{3}','27 ^{\circ}C, 200 \mug/m^{3}','27 ^{\circ}C, 250 \mug/m^{3}');
xlim([0 80]);


% Probability of match lasting three sets
% SELECT: winProbDiff_dis, HighCost_Ability_dis OR cl_ch_ratio along the horizontal axis
figure;
subplot(1,1,1);
plot(winProbDiff_dis,prob3setter(:,1),'g',winProbDiff_dis,prob3setter(:,2),'--r',winProbDiff_dis,prob3setter(:,3),'--ro',winProbDiff_dis,prob3setter(:,4),':k',winProbDiff_dis,prob3setter(:,5),':ko','LineWidth',3);
%plot(HighCost_Ability_dis,prob3setter(:,1),'g',HighCost_Ability_dis,prob3setter(:,2),'--r',HighCost_Ability_dis,prob3setter(:,3),'--ro',HighCost_Ability_dis,prob3setter(:,4),':k',HighCost_Ability_dis,prob3setter(:,5),':ko','LineWidth',3);
%plot(cl_ch_ratio    ,prob3setter(:,1),'g',cl_ch_ratio    ,prob3setter(:,2),'--r',cl_ch_ratio    ,prob3setter(:,3),'--ro',cl_ch_ratio    ,prob3setter(:,4),':k',cl_ch_ratio    ,prob3setter(:,5),':ko','LineWidth',3);
ylabel('Probability of match lasting three battles');
xlabel('Difference in opponents'' pre-match winning probability (%)');
%xlabel('High-cost player''s rank (competing with player ranked 1)');
%xlabel('Players'' relative marginal costs, c_l/c_h');
legend('27 ^{\circ}C, 150 \mug/m^{3}','32 ^{\circ}C, 150 \mug/m^{3}','37 ^{\circ}C, 150 \mug/m^{3}','27 ^{\circ}C, 200 \mug/m^{3}','27 ^{\circ}C, 250 \mug/m^{3}');
xlim([0 80]);


% Match outcome predictability as environmental quality varies
% SELECT: winProbDiff_dis, HighCost_Ability_dis OR cl_ch_ratio along the horizontal axis
figure;
subplot(1,1,1);
plot(winProbDiff_dis,winProb_l_dis(:,1),'g',winProbDiff_dis,winProb_l_dis(:,2),'--r',winProbDiff_dis,winProb_l_dis(:,3),'--ro',winProbDiff_dis,winProb_l_dis(:,4),':k',winProbDiff_dis,winProb_l_dis(:,5),':ko','LineWidth',3);
%plot(HighCost_Ability_dis,winProb_l_dis(:,1),'g',HighCost_Ability_dis,winProb_l_dis(:,2),'--r',HighCost_Ability_dis,winProb_l_dis(:,3),'--ro',HighCost_Ability_dis,winProb_l_dis(:,4),':k',HighCost_Ability_dis,winProb_l_dis(:,5),':ko','LineWidth',3);
%plot(cl_ch_ratio    ,winProb_l_dis(:,1),'g',cl_ch_ratio    ,winProb_l_dis(:,2),'--r',cl_ch_ratio    ,winProb_l_dis(:,3),'--ro',cl_ch_ratio    ,winProb_l_dis(:,4),':k',cl_ch_ratio    ,winProb_l_dis(:,5),':ko','LineWidth',3);
ylabel('Low-cost player''s match winning probability');
xlabel('Difference in opponents'' pre-match winning probability (%)');
%xlabel('High-cost player''s rank (competing with player ranked 1)');
%xlabel('Players'' relative marginal costs, c_l/c_h');
legend('27 ^{\circ}C, 150 \mug/m^{3}','32 ^{\circ}C, 150 \mug/m^{3}','37 ^{\circ}C, 150 \mug/m^{3}','27 ^{\circ}C, 200 \mug/m^{3}','27 ^{\circ}C, 250 \mug/m^{3}');
xlim([0 80]);
ylim([.5 1]); 


% Players' expected payoffs as environmental quality varies; EDIT plot command & legend if delta_no_gridpoints changes
% SELECT: winProbDiff_dis, HighCost_Ability_dis OR cl_ch_ratio along the horizontal axis
figure;
subplot(1,1,1);
plot(winProbDiff_dis,exppayoff_l(:,1),'g',winProbDiff_dis,exppayoff_l(:,2),'--r',winProbDiff_dis,exppayoff_l(:,3),'--ro',winProbDiff_dis,exppayoff_l(:,4),':k',winProbDiff_dis,exppayoff_l(:,5),':ko','LineWidth',3);
%plot(HighCost_Ability_dis,exppayoff_l(:,1),'g',HighCost_Ability_dis,exppayoff_l(:,2),'--r',HighCost_Ability_dis,exppayoff_l(:,3),'--ro',HighCost_Ability_dis,exppayoff_l(:,4),':k',HighCost_Ability_dis,exppayoff_l(:,5),':ko','LineWidth',3);
%plot(cl_ch_ratio    ,exppayoff_l(:,1),'g',cl_ch_ratio    ,exppayoff_l(:,2),'--r',cl_ch_ratio    ,exppayoff_l(:,3),'--ro',cl_ch_ratio    ,exppayoff_l(:,4),':k',cl_ch_ratio    ,exppayoff_l(:,5),':ko','LineWidth',3);
ylabel('Low-cost player''s expected payoff (1000 $)');
xlabel('Difference in opponents'' pre-match winning probability (%)');
%xlabel('High-cost player''s rank (competing with player ranked 1)');
%xlabel('Players'' relative marginal costs, c_l/c_h');
legend('27 ^{\circ}C, 150 \mug/m^{3}','32 ^{\circ}C, 150 \mug/m^{3}','37 ^{\circ}C, 150 \mug/m^{3}','27 ^{\circ}C, 200 \mug/m^{3}','27 ^{\circ}C, 250 \mug/m^{3}');
xlim([0 80]);
%ylim([0 12]);

figure;
subplot(1,1,1);
plot(winProbDiff_dis,exppayoff_h(:,1),'g',winProbDiff_dis,exppayoff_h(:,2),'--r',winProbDiff_dis,exppayoff_h(:,3),'--ro',winProbDiff_dis,exppayoff_h(:,4),':k',winProbDiff_dis,exppayoff_h(:,5),':ko','LineWidth',3);
%plot(HighCost_Ability_dis,exppayoff_h(:,1),'g',HighCost_Ability_dis,exppayoff_h(:,2),'--r',HighCost_Ability_dis,exppayoff_h(:,3),'--ro',HighCost_Ability_dis,exppayoff_h(:,4),':k',HighCost_Ability_dis,exppayoff_h(:,5),':ko','LineWidth',3);
%plot(cl_ch_ratio    ,exppayoff_h(:,1),'g',cl_ch_ratio    ,exppayoff_h(:,2),'--r',cl_ch_ratio    ,exppayoff_h(:,3),'--ro',cl_ch_ratio    ,exppayoff_h(:,4),':k',cl_ch_ratio    ,exppayoff_h(:,5),':ko','LineWidth',3);
ylabel('High-cost player''s expected payoff (1000 $)');
xlabel('Difference in opponents'' pre-match winning probability (%)');
%xlabel('High-cost player''s rank (competing with player ranked 1)');
%xlabel('Players'' relative marginal costs, c_l/c_h');
legend('27 ^{\circ}C, 150 \mug/m^{3}','32 ^{\circ}C, 150 \mug/m^{3}','37 ^{\circ}C, 150 \mug/m^{3}','27 ^{\circ}C, 200 \mug/m^{3}','27 ^{\circ}C, 250 \mug/m^{3}');
xlim([0 80]);

%}












%{
% Plot Figure 6 on best response functions for dynamic contest model section (no estimates here). k is fixed at 1. Symmetric win-loss prize spread. 
BR_no_gridpoints=120;
xrival=zeros(BR_no_gridpoints,1);
BR_c_point5 =zeros(BR_no_gridpoints,1);
BR_c_point25=zeros(BR_no_gridpoints,1);
for i=1:BR_no_gridpoints;
    xrival(i)=.01+(i-1)*.01;
    BR_c_point5(i) =(xrival(i)/.5 )^.5-xrival(i);
    BR_c_point25(i)=(xrival(i)/.25)^.5-xrival(i);
end;
    %Add rival effort closer to 0
    xrival=[.000001; .0001; xrival];
    BR_c_point5 =[(.000001/.5 )^.5-.000001; (.0001/.5 )^.5-.0001; BR_c_point5];
    BR_c_point25=[(.000001/.25)^.5-.000001; (.0001/.25)^.5-.0001; BR_c_point25];
figure;
subplot(1,2,1);
plot(BR_c_point5,xrival,':r',xrival,BR_c_point5,'-r*','LineWidth',5);
title('Optimal effort levels for symmetric case')
xlabel('Low-cost player''s effort, x_l');
ylabel('High-cost player''s effort, x_h');
legend('Low-cost player''s best response to rival''s choice','High-cost player''s best response to rival''s choice');
subplot(1,2,2);
plot(BR_c_point25,xrival,':b',xrival,BR_c_point5,'-r*','LineWidth',5);
title('Optimal effort levels for asymmetric case')
xlabel('Low-cost player''s effort, x_l');
ylabel('High-cost player''s effort, x_h');
legend('Low-cost player''s best response to rival''s choice','High-cost player''s best response to rival''s choice');
%}



%{
% Plot Figures 7, A.1, A.3, and A.4
no_gridpoints=100;
k=1;
%k=.7;    % Change randomness here   
ones=ones(no_gridpoints,1);
[p_llose3_symm, p_llose2_lwon1_symm, p_llose2_llost1_symm, p_llose1_symm, xl3_symm, xh3_symm, xl2_lwon1_symm, xh2_lwon1_symm, xl2_llost1_symm, xh2_llost1_symm, xl1_symm, xh1_symm, delta_plot] = plot_transitions_effort_over_delta(no_gridpoints,.5,.5,1,1,k);     % Inputs are no_gridpoints, c_l, c_h, V_l, V_h, k   
[p_llose3_asymm,p_llose2_lwon1_asymm,p_llose2_llost1_asymm,p_llose1_asymm,xl3_asymm,xh3_asymm,xl2_lwon1_asymm,xh2_lwon1_asymm,xl2_llost1_asymm,xh2_llost1_asymm,xl1_asymm,xh1_asymm,delta_plot] = plot_transitions_effort_over_delta(no_gridpoints,.25,.5,1,1,k);
%[p_llose3_asymm,p_llose2_lwon1_asymm,p_llose2_llost1_asymm,p_llose1_asymm,xl3_asymm,xh3_asymm,xl2_lwon1_asymm,xh2_lwon1_asymm,xl2_llost1_asymm,xh2_llost1_asymm,xl1_asymm,xh1_asymm,delta_plot] = plot_transitions_effort_over_delta(no_gridpoints,.5,1,1,1,k);
%[p_llose3_veryasymm,p_llose2_lwon1_veryasymm,p_llose2_llost1_veryasymm,p_llose1_veryasymm,xl3_veryasymm,xh3_veryasymm,xl2_lwon1_veryasymm,xh2_lwon1_veryasymm,xl2_llost1_veryasymm,xh2_llost1_veryasymm,xl1_veryasymm,xh1_veryasymm,delta_plot] = plot_transitions_effort_over_delta(no_gridpoints,.25,1,1,1,k);

% Figure 7 on transitions for dynamic contest model section (no estimates here). k is fixed at 1. Symmetric contest prizes.   
figure;
subplot(1,2,1);
plot(delta_plot,p_llose1_symm,'r',delta_plot,p_llose2_llost1_symm,':r',delta_plot,p_llose2_lwon1_symm,'-r*','LineWidth',5);
title('Stage transition probabilities for symmetric case')
xlabel('Environmental disutility parameter (per battle)');
ylabel('Transition probabilities');
legend('Player wins battle 1 or wins battle 3','Player wins battle 2 conditional on battle-1 win','Player wins battle 2 conditional on battle-1 loss');
ylim([0 1]);
subplot(1,2,2);
plot(delta_plot,p_llose1_asymm,'r',delta_plot,p_llose2_llost1_asymm,':r',delta_plot,p_llose2_lwon1_asymm,'-r*',delta_plot,p_llose3_asymm,'--dr','LineWidth',5);
title('Stage transition probabilities for asymmetric case')
xlabel('Environmental disutility parameter (per battle)');
ylabel('Transition probabilities');
legend('High-cost player wins battle 1','High-cost player wins battle 2 | battle-1 win','High-cost player wins battle 2 | battle-1 loss','High-cost player wins battle 3');
ylim([0 1]);

% Figure A.1 on contest winning probabilities for dynamic contest model section (no estimates here). Symmetric contest prizes.  
% Sum of probabilities of ll, lhl, hll (corresponding to the three additive terms in this order). Note winProb_l_symm=winProb_h_symm
winProb_l_symm  = (ones-p_llose1_symm) .*(ones-p_llose2_lwon1_symm)  + (ones-p_llose1_symm) .*p_llose2_lwon1_symm  .*(ones-p_llose3_symm)  + p_llose1_symm .*(ones-p_llose2_llost1_symm) .*(ones-p_llose3_symm);
winProb_l_asymm = (ones-p_llose1_asymm).*(ones-p_llose2_lwon1_asymm) + (ones-p_llose1_asymm).* p_llose2_lwon1_asymm.*(ones-p_llose3_asymm) + p_llose1_asymm.*(ones-p_llose2_llost1_asymm).*(ones-p_llose3_asymm);
% Sum of probabilities of lhh, hlh, hh (corresponding to the three additive terms in this order)
winProb_h_asymm = (ones-p_llose1_asymm).*p_llose2_lwon1_asymm.*p_llose3_asymm + p_llose1_asymm.*(ones-p_llose2_llost1_asymm).*p_llose3_asymm + p_llose1_asymm.*p_llose2_llost1_asymm;
figure;
subplot(1,2,1);
plot(delta_plot,winProb_l_symm,'r','LineWidth',5);
title('Contest winning probability for symmetric case')
xlabel('Environmental disutility parameter (per battle)');
ylabel('Contest winning probability');
legend('Player''s contest winning probability');
ylim([0 .9]);
subplot(1,2,2);
plot(delta_plot,winProb_l_asymm,':b',delta_plot,winProb_h_asymm,'r','LineWidth',5);
title('Contest winning probability for asymmetric case')
xlabel('Environmental disutility parameter (per battle)');
ylabel('Contest winning probability');
legend('Low-cost player''s contest winning probability','High-cost player''s contest winning probability');
ylim([0 .9]);

% Figure A.3 on model-predicted probability of reaching a third set for dynamic contest model section (no estimates here). Symmetric contest prizes. 
% Sum of probabilities of lhl, lhh, hll, hlh (corresponding to the three additive terms in this order)
probThreeSetter_symm  = (ones-p_llose1_symm) .*p_llose2_lwon1_symm .*(ones-p_llose3_symm)  + (ones-p_llose1_symm) .*p_llose2_lwon1_symm .*p_llose3_symm  + p_llose1_symm .*(ones-p_llose2_llost1_symm) .*(ones-p_llose3_symm)  + p_llose1_symm .*(ones-p_llose2_llost1_symm) .*p_llose3_symm;
probThreeSetter_asymm = (ones-p_llose1_asymm).*p_llose2_lwon1_asymm.*(ones-p_llose3_asymm) + (ones-p_llose1_asymm).*p_llose2_lwon1_asymm.*p_llose3_asymm + p_llose1_asymm.*(ones-p_llose2_llost1_asymm).*(ones-p_llose3_asymm) + p_llose1_asymm.*(ones-p_llose2_llost1_asymm).*p_llose3_asymm;
figure;
subplot(1,2,1);
plot(delta_plot,probThreeSetter_symm,'r','LineWidth',5);
title('Probability of 3-battle contest for symmetric case')
xlabel('Environmental disutility parameter (per battle)');
ylabel('Probability of 3-battle contest');
legend('Probability of 3-battle contest');
ylim([0 .4]);
subplot(1,2,2);
plot(delta_plot,probThreeSetter_asymm,'r','LineWidth',5);
title('Probability of 3-battle contest for asymmetric case')
xlabel('Environmental disutility parameter (per battle)');
ylabel('Probability of 3-battle contest');
legend('Probability of 3-battle contest');
ylim([0 .4]);

% Figure A.4 on effort levels in stage 1 for dynamic contest model section (no estimates here). k is fixed at 0.7. Symmetric contest prizes.  
figure;
subplot(1,2,1);
plot(delta_plot,xl1_symm,'r','LineWidth',5);
title('Battle-1 effort levels for symmetric case')
xlabel('Environmental disutility parameter (per battle)');
ylabel('Optimal effort');
legend('Player''s effort');
ylim([0 .5]);
subplot(1,2,2);
plot(delta_plot,xl1_asymm,':b',delta_plot,xh1_asymm,'r','LineWidth',5);
title('Battle-1 effort levels for asymmetric case')
xlabel('Environmental disutility parameter (per battle)');
ylabel('Optimal effort');
legend('Low-cost player''s effort','High-cost player''s effort');
ylim([0 .5]);

% Figure A.4 on effort levels in stage 2 for dynamic contest model section (no estimates here). k is fixed at 0.7. Symmetric contest prizes.
figure;
subplot(1,2,1);
plot(delta_plot,xl2_lwon1_symm,'r',delta_plot,xl2_llost1_symm,'-r*','LineWidth',5);
title('Battle-2 effort levels for symmetric case')
xlabel('Environmental disutility parameter (per battle)');
ylabel('Optimal effort');
legend('Player''s effort | battle-1 win','Player''s effort | battle-1 loss');
ylim([0 .35]);
subplot(1,2,2);
plot(delta_plot,xl2_lwon1_asymm,':b',delta_plot,xl2_llost1_asymm,'-b*',delta_plot,xh2_llost1_asymm,'r',delta_plot,xh2_lwon1_asymm,'-r*','LineWidth',5);
title('Battle-2 effort levels for asymmetric case')
xlabel('Environmental disutility parameter (per battle)');
ylabel('Optimal effort');
legend('Low-cost player''s effort | battle-1 win','Low-cost player''s effort | battle-1 loss','High-cost player''s effort | battle-1 win','High-cost player''s effort | battle-1 loss');
ylim([0 .35]);

% Figure A.4 on effort levels in stage 3 for dynamic contest model section (no estimates here). k is fixed at 0.7. Symmetric contest prizes.
figure;
subplot(1,2,1);
plot(delta_plot,xl3_symm,'r','LineWidth',5);
title('Battle-3 effort levels for symmetric case')
xlabel('Environmental disutility parameter (per battle)');
ylabel('Optimal effort');
legend('Player''s effort');
ylim([0 .7]);
subplot(1,2,2);
plot(delta_plot,xl3_asymm,':b',delta_plot,xh3_asymm,'r','LineWidth',5);
title('Battle-3 effort levels for asymmetric case')
xlabel('Environmental disutility parameter (per battle)');
ylabel('Optimal effort');
legend('Low-cost player''s effort','High-cost player''s effort');
ylim([0 .7]);
%}



%{
% Plot Figure A.2 on match winning probability for stronger player against cl_ch_ratio ranging from asymmetry at .2 to symmetry at 1 (extremities for 81 gridpoints). Symmetric contest prizes.
cost_no_gridpoints=81;      %from 1 to 0.2 in steps of 0.01
%Plot for delta=0, delta=.05, and delta=.1 (elements 1, 51 & 101 in vector with 101 gridpoints)
delta_no_gridpoints=101;
k=1;
%k=.7;    % Change randomness here
cl_ch_ratio       =zeros(cost_no_gridpoints,1);
winProb_l_env_good=zeros(cost_no_gridpoints,1);
winProb_l_env_med =zeros(cost_no_gridpoints,1);
winProb_l_env_poor=zeros(cost_no_gridpoints,1);

for i=1:cost_no_gridpoints;
    cl_ch_ratio(i)=(100-cost_no_gridpoints+1)/100+(i-1)*.01;
    [p_llose3,p_llose2_lwon1,p_llose2_llost1,p_llose1,~,~,~,~,~,~,~,~,delta] = plot_transitions_effort_over_delta(delta_no_gridpoints,cl_ch_ratio(i),1,1,1,k);     % Inputs are no_gridpoints, c_l, c_h, V_l, V_h, k
    winProb_l_env_good(i)=(1-p_llose1(1))  *(1-p_llose2_lwon1(1))   + (1-p_llose1(1))  *p_llose2_lwon1(1)  *(1-p_llose3(1))   + p_llose1(1)  *(1-p_llose2_llost1(1))  *(1-p_llose3(1));
    winProb_l_env_med(i) =(1-p_llose1(51)) *(1-p_llose2_lwon1(51))  + (1-p_llose1(51)) *p_llose2_lwon1(51) *(1-p_llose3(51))  + p_llose1(51) *(1-p_llose2_llost1(51)) *(1-p_llose3(51));
    winProb_l_env_poor(i)=(1-p_llose1(end))*(1-p_llose2_lwon1(end)) + (1-p_llose1(end))*p_llose2_lwon1(end)*(1-p_llose3(end)) + p_llose1(end)*(1-p_llose2_llost1(end))*(1-p_llose3(end));
end;

figure;
subplot(1,1,1);
plot(cl_ch_ratio,winProb_l_env_good,'k',cl_ch_ratio,winProb_l_env_med,':k',cl_ch_ratio,winProb_l_env_poor,'-k*','LineWidth',5);
title('Contest outcome predictability as environmental quality shifts')
xlabel('Players'' relative marginal costs, c_l/c_h');
ylabel('Low-cost player''s contest winning probability');
legend('High environmental quality, \delta = 0','Medium environmental quality, \delta = 0.05','Low environmental quality, \delta = 0.1');
ylim([0.5 1]);
%}


